*put path and dofile in the same folder
*************************set path and use QOG data********************
cd "path_of_dofile_and_datafile"
use qog_patent_vdem_merged,clear
*************************data cleaning***************************
gen hl20= wdi_incsh20h/wdi_incsh20l
label var hl20 "The ratio of income share held by highest 20% to lowest 20% " 
gen middle=wdi_incsh203+wdi_incsh204
gen l_patent=log(patent)
label var van_part "Political participation"
label var wdi_unempilo "Unemployment, total (% of total labor force)"
label var wdi_gini "Gini index (World Bank estimate)"
label var  wdi_internet "Internet penetration"
label var middle "middle class defined by 3rd and 4th 20% income group"
label var l_patent "The log of patent grants"
//north-south-oecd
gen lower=wdi_incsh203+wdi_incsh204+wdi_incsh20l
label var north "global north counries"
gen oecd=1 if oecd_sizegdp_t1!=.
gen p46=(wdi_incsh20h+wdi_incsh202)/(wdi_incsh204+wdi_incsh20l)
gen p28=wdi_incsh20h/wdi_incsh20l
replace oecd=0 if oecd!=1
label var oecd "oecd countries"
/*keep wbgi_cce wdi_gerp wdi_gert $y wdi_internet van_part   ///
 wdi_unempilo  wdi_empind wdi_poprulgr wdi_agedr ccode year ccodealp ///
wdi_incsh20h wdi_incsh20l hl20 wdi_mobile patent wdi_gdpcapgr oecd north l_patent wdi_unempedui ///
wdi_incsh203 wdi_incsh204 middle van_index lower*/
//data processing
gen wdi_2=wdi_internet^2

global y p46
replace wdi_pop=wdi_pop/1000000
global xlist ictd_taxinsc wdi_gdpcapcur wdi_agedr wdi_popgr wdi_pop wdi_gdpind wdi_gdpagr wdi_trade
global mlist wbgi_cce wdi_gerp wdi_povgap365 wdi_unempyne 
global rlist wdi_incsh20h wdi_incsh20l wdi_gini wdi_mobile patent
xtreg $y wdi_internet $xlist
gen a=e(sample)
bysort a:gen count=_N
drop if count<4 //drop the countries with fewer observation
save ict_ineaquality,replace

************************main text analysis (Table 2)***************

use ict_ineaquality,clear
xtset ccode year
//data description
sum $y wdi_internet $xlist $mlist $rlist
des $y wdi_internet $xlist $mlist $rlist

xtreg $y wdi_mobile i.year if year>=2017,fe 

//basic results 
reg $y wdi_internet $xlist 
est store reg1
xtgls $y wdi_internet $xlist i.year,panels(hetero)
est store reg2
xtreg $y wdi_internet i.year,re 
est store reg3
xtreg $y wdi_internet $xlist i.year,fe 
est store reg4
xtreg $y wdi_internet wdi_2 $xlist i.year ,fe 
est store reg5
xtreg $y L(1/2).$y wdi_internet wdi_2 $xlist v2mecenefi v2smregcon i.year, fe 
est store reg6
esttab  reg1 reg2 reg3 reg4 reg5 reg6 , cells(b(star fmt(4)) t(par fmt(4))) /// 
legend label varwidth(60) varlabels(_cons constant) mtitles("FE" "FE") replace /// 
stats(N r2 N_g , fmt(0 4) label(Observation R-squared )) star(* 0.1 ** 0.05 *** 0.01) /*table 3*/
capture hausman reg3 reg4,sigmamore
esttab  reg1 reg2 reg3 reg4 reg5 reg6 using baseline.rtf, cells(b(star fmt(4)) t(par fmt(4))) /// 
legend label varwidth(60) varlabels(_cons constant) mtitles("FE" "FE") replace /// 
stats(N r2 N_g , fmt(0 4) label(Observation R-squared )) star(* 0.1 ** 0.05 *** 0.01) /*table 3*/
capture hausman reg3 reg4,sigmamore

****************************moderating effect (Table 3)
des wbgi_cce wdi_gerp wdi_gert

//
sum wbgi_cce,de
gen level1=.
replace level1=0 if wbgi_cce<=r(mean)
replace level1=1 if wbgi_cce>r(mean)&wbgi_cce<.
gen cce_inter=level1*wdi_internet
xtreg $y wdi_internet $xlist i.year if level1==1,fe
est store reg1
xtreg $y wdi_internet $xlist i.year if level1==0,fe
est store reg2
xtreg wdi_gini cce_inter level1 wdi_internet $xlist i.year,fe
est store reg3

//
sum wdi_gerp,de
gen level2=.
replace level2=0 if wdi_gerp<=r(mean)
replace level2=1 if wdi_gerp>r(mean)&wdi_gerp<.
gen gerp_inter=level2*wdi_internet
xtreg $y wdi_internet  $xlist i.year if level2==1,fe
est store reg4
xtreg $y wdi_internet  $xlist i.year if level2==0,fe
est store reg5
xtreg $y wdi_internet gerp_inter level2  $xlist i.year,fe
est store reg6
esttab reg1 reg2 reg3 reg4 reg5 reg6, cells(b(star fmt(4)) t(par fmt(4))) /// 
legend label varwidth(60) varlabels(_cons Constant) mtitles("High" "Low" "Interaction" "High" "Low" "Interaction") /// 
keep (wdi_internet  gerp_inter cce_inter  level2  level1   _cons) ///
stats(N r2 N_g _w, fmt(0 4) label(Observation R-squared )) star(* 0.1 ** 0.05 *** 0.01) /*table 3*/

esttab reg1 reg2 reg3 reg4 reg5 reg6 using moderation.rtf, cells(b(star fmt(4)) t(par fmt(4))) /// 
legend label varwidth(60) varlabels(_cons Constant) mtitles("High" "Low" "Interaction" "High" "Low" "Interaction") /// 
keep (wdi_internet  gerp_inter cce_inter  level2  level1   _cons) replace ///
stats(N r2 N_g _w, fmt(0 4) label(Observation R-squared )) star(* 0.1 ** 0.05 *** 0.01) /*table 3*/


************************************mediating effect(Table 4)
des wdi_empind wdi_gini wdi_gini wdi_gini
est store reg
winsor wdi_povgap365, p(0.1) gen(wpovgap)
winsor wdi_unempedui, p(0.05) gen(wunemp)
winsor wdi_empind,p(0.1) gen(wempind)
gen lnwdi_povgap365=log(wdi_povgap365)
gen lnempind=log(wdi_empind)
gen lnempedui=log(wdi_unempedui)
xtreg wpovgap wdi_internet $xlist i.year,fe
est store reg1
xtreg $y wpovgap $xlist i.year,fe
est store reg2
xtreg wdi_unempyne wdi_internet $xlist i.year,fe
est store reg3
xtreg $y wdi_unempyne $xlist i.year,fe
est store reg4
esttab  reg1 reg2 reg3 reg4 using mediating.rtf,replace  cells(b(star fmt(4)) t(par fmt(4)))  /// 
legend label varwidth(60) varlabels(_cons constant) mtitles("FE" "FE") /// 
stats(N r2 N_g , fmt(0 4) label(Observation R-squared )) star(* 0.1 ** 0.05 *** 0.01) /*table 3*/

************************************robust test(Table 5)

xi i.year
forv i = 2000/2019{
	rename _Iyear_`i' yeard`i'
}
drop yeard2000
xtreg wdi_incsh20h wdi_internet $xlist i.year,fe 
est store reg1
xtreg wdi_incsh20l wdi_internet $xlist i.year,fe 
est store reg2
xtreg wdi_gini wdi_internet $xlist i.year,fe 
est store reg3
 xtivreg $y (wdi_internet=wdi_mobile l_patent) $xlist yeard*,fe first
est store reg4
xtoverid
dmexogxt
esttab  reg1 reg2 reg3 reg4 using robust.rtf,  cells(b(star fmt(4)) t(par fmt(4)))  /// 
legend label varwidth(60) varlabels(_cons constant) mtitles("FE" "FE" ) replace /// 
stats(N r2 N_g , fmt(0 4) label(Observation R-squared )) star(* 0.1 ** 0.05 *** 0.01) /*table 3*/


************************************spatial vairation (Figure 3)
global size size(*1.2)
global scale scale(1)
xtreg $y c.wdi_internet wdi_2 $xlist i.year ,fe
est store all
 . global tp: dis %2.0f -_b[wdi_internet]/(2*_b[wdi_2]) //turning point
 . dis "Turn point = " $tp 
 . local b0 = _b[_cons]
 . local b1 = _b[wdi_internet]
 . local b2 = _b[wdi_2]
 . sum wdi_internet
 . local min:  dis %3.0f r(min)
 . local max:  dis %3.0f r(max)
 . twoway ( function y = `b2'*x^2 + `b1'*x + `b0', ///
           range(0 100) ), ///
  	 ytitle("Income inequality") xtitle("Digitalization") ///
  	 xline($tp,  lp(dash) lc(red)) ///
  	 xlabel($tp "($tp)" 0(50)100) $scale ///
  	 legend(off) scheme(s1mono) saving(all,replace) ///
	 b1title("(a) all countries", `size') 
  . graph export "all.png", replace  //保存图形

xtreg $y c.wdi_internet wdi_2 $xlist i.year if euro==1,fe 
est store origin
 . global tp: dis %2.0f -_b[wdi_internet]/(2*_b[wdi_2]) //turning point
 . dis "Turn point = " $tp 
 . local b0 = _b[_cons]
 . local b1 = _b[wdi_internet]
 . local b2 = _b[wdi_2]
 . sum wdi_internet
 . local min:  dis %3.0f r(min)
 . local max:  dis %3.0f r(max)
 . twoway ( function y = `b2'*x^2 + `b1'*x + `b0', ///
           range(0 100) ), ///
  	 ytitle("Income inequality") xtitle("Digitalization") ///
  	 xline($tp,  lp(dash) lc(red)) ///
  	 xlabel($tp "($tp)" 0(50)100) $scale ///
  	 legend(off) scheme(s1mono) saving(origin,replace) ///
	 b1title("(b) origin of digitalization", `size') 
  . graph export "Origin.png", replace  //保存图形

  xtreg $y c.wdi_internet wdi_2 $xlist i.year if euro==0,fe 
  est store influenced
  . global tp: dis %2.0f -_b[wdi_internet]/(2*_b[wdi_2]) //turning point
 . dis "Turn point = " $tp 
 . local b0 = _b[_cons]
 . local b1 = _b[wdi_internet]
 . local b2 = _b[wdi_2]
 . sum wdi_internet
 . local min:  dis %3.0f r(min)
 . local max:  dis %3.0f r(max)
 . twoway ( function y = `b2'*x^2 + `b1'*x + `b0', ///
           range(0 100) ), ///
  	 ytitle("Income inequality") xtitle("Digitalization") ///
  	 xline($tp,  lp(dash) lc(red)) ///
  	 xlabel($tp "($tp)" 0(50)100) $scale ///
  	 legend(off) scheme(s1mono) saving(influenced,replace) ///
	 b1title("(c) influenced by digitalization", `size') 
  . graph export "Influenced.png", replace  
  
xtreg $y c.wdi_internet wdi_2 $xlist i.year if north==0,fe 
est store south
 . global tp: dis %2.0f -_b[wdi_internet]/(2*_b[wdi_2]) //turning point
 . dis "Turn point = " $tp 
 . local b0 = _b[_cons]
 . local b1 = _b[wdi_internet]
 . local b2 = _b[wdi_2]
 . sum wdi_internet
 . local min:  dis %3.0f r(min)
 . local max:  dis %3.0f r(max)
 . twoway ( function y = `b2'*x^2 + `b1'*x + `b0', ///
           range(0 100) ), ///
  	 ytitle("Income inequality") xtitle("Digitalization") ///
  	 xline($tp,  lp(dash) lc(red)) ///
  	 xlabel($tp "($tp)" 0(50)100) $scale ///
  	 legend(off) scheme(s1mono) saving(global_south,replace)  ///
	 b1title("(d) global_south countries", `size') 
  . graph export "south.png", replace  

graph combine all.gph origin.gph influenced.gph global_south.gph
graph export combine.emf, replace 
graph export combine.png, replace 

esttab all origin influenced south using spatial_variance.rtf,cells(b(star fmt(4)) t(par fmt(4)))/// 
legend label varwidth(60) varlabels(_cons constant) mtitles("all" "origin" "influenced" "south") /// 
 replace ///
stats(N N_g r2, fmt(0 0 4) label(Observation Countries R-squared )) star(* 0.1 ** 0.05 *** 0.01) 
